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ABSTRACT 


The objective of this study is to predict the pressure 
response of a saturated liquid-vapor system when undergoing a 

venting or depressurization process in zero gravity at low vent 
rates. An experimental investigation of the venting of cylindri- 
cal containers partially filled with initially saturated liquids 
was previously conducted under zero-gravity conditions at the 
NASA Lewis Research Center 5-second zero gravity facility, and 
compared with an analytical model which incorporated the effect 
of interfacial mass transfer on the ullage pressure response 
during venting. A new model is presented here to improve the 
estimation of the interfacial mass transfer. Duhammel's super- 
position integral is incorporated to approximate the transient 
temperature response of the interface, treating the liquid as a 
semi-infinite solid with conduction heat transfer. Account is 
also taken of the condensation taking place within the bulk of a 
saturated vapor as isentropic expansion takes place. Computa- 
tional results are presented for the venting of R-ll from a given 
vessel and initial state for five different venting rates over a 
period of three seconds, and are compared to the prior NASA ex- 
periments. An improvement in the prediction of the final pres- 
sure takes place, but is still considerably below the measure- 
ments. This is attributed to neglecting the evaporation taking 
place at the meniscus . 



I . INTRODUCTION 


This is a report of an analytic study of the venting pro- 
cess, under microgravity, of a container with liquid-vapor 

phases. Pure substances only are considered here. It may be 
desireable to extend the analysis later to mixtures of substances 
having different boiling points should the possibility arise for 
the future storage of such mixtures in space. 

Microgravity indicates the drastic reduction of body forces 
which, in turn, implies the drastic reduction of natural convec- 
tion motion associated with temperature differences within the 
fluid(s). With liquid-vapor phases present, this also implies 
that effects associated with surface tension may become signifi- 
cant. These effects include the absence of a flat interface, the 
possible presence of thermophorysis , and the variable location of 
the vapor volume. This latter is important for the venting 
process in practice, since the process desired must be specified, 
whether venting of pure vapor, pure liquid, or a mixture of both. 
For purposes of the present study surface tension is considered 
only insofar as it affects the amount of liquid-vapor interfacial 
area and the contact angle made at the solid-liquid-vapor contact 
region. The curvature may be neglected in the description of the 
temperature distribution in the liquid near the liquid-vapor 
interface, except where very small vapor bubbles are present. 

Such would mean that vapor nucleation and boiling are taking 
place , which are also excluded from consideration at present . 

The objective of the study is to model the venting process 
so as to be able to predict the pressure-time behavior within the 
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container , incorporating the relevent physical mechanisms . In 
the broadest sense the process is similar to the problem of de- 
scribing the blowdown taking place in a nuclear power plant loss- 
of-coolant accident (LOCA). Because of its importance to the 
safety of nuclear power plants many studies of this have been 
conducted, and an abundance of literature exists. However, very 
little is directly applicable to the present problem. 

The LOCA is generally an uncontrolled process in that the 
circumstance as to the location and size of the discharge 
opening is not known in advance. The fluid discharged may be 
liquid only, vapor only, or a mixture of the two, depending on 
the location and size of the opening. If the opening is suf- 
ficiently small the pressure decrease may be described reasonably 
well by assuming saturation properties within the vessel, pro- 
vided the two-phase critical flow occuring at the opening is 
properly modelled. Some difficulties still exist in achieving 
this [1,2] without introducing some degree of empiricism for 
the discharge coefficient. If the opening is relatively large 
the pressure drops so rapidly that bulk vapor nucleation occurs 
with a two-phase discharge. The upper limit of this is the pipe 
blowdown problem: A one-dimensional analysis and confirming 

experiments of this process, beginning with a subcooled liquid, 
were conducted [3] in which the time for complete discharge was 
on the order of 200-300 ms. 

The zero gravity venting problem under consideration here , 

on the other hand, is expected to be a more orderly and planned 
process . Nevertheless , the description of the physical process 
can become relatively complex when the transient behavior in the 
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vicinity of the liquid-vapor interface is taken into considera- 
tion . An analysis of the adiabatic venting of the homogeneous 
two-phase contents of a vessel was conducted, assuming equili- 
brium saturation conditions [4], A special model for choked 
two-phase flow at the discharge opening was derived and used. 
Since equilibrium conditions were assumed, the size and location 
of the vapor during the process was irrelevant. In a later work 
[5), the same author attempted to include vapor generation asso- 
ciated with external heat transfer in the process, retaining the 
assumption of internal thermodynamic equilibrium. In addition, 
various venting locations were considered: Vapor venting only 

(top vapor blowdown), liquid venting only (bottom liquid blow- 
down), and bottom mixture blowdown. 

The zero gravity venting problem at hand is concerned only 
with the venting of vapor, either with a known vent rate or vent 
size. It is further considered that the venting rate is suffi- 
ciently slow that no vapor nucleating sites are formed, so that 
any phase changes taking place occur only at existing interfaces. 
This places an upper limit on the rate of pressure drop permis- 
sible [6,7,8]. 

It has been demonstrated that models assuming thermodynamic 
equilibrium can give reasonable results only with the inclusion 
of sufficient empirical coefficients, since non-equilibrium con- 
ditions in fact are present. The challenge is to minimize the 
number of empirical coefficients, using non-equilibrium analysis 
only for those aspects of the process where it has a significant 
role . 


- 4 - 



The work reported below is an extension of the analysis 

presented in [ 9 ] , employing a different model for the estimation 

of the mass transfer at the liquid-vapor interface and in includ- 
ing the phase change associated with the expansion of the satur- 
ated vapor as the venting process takes place. Discrepancies 
still exist between the computations conducted here and the ex- 
perimental results reported in [9], with the model predicting a 
greater pressure decrease than is measured. The transient temp- 
erature at the liquid-vapor interface is computed by incorporat- 
ing Duhammel's superposition integral in the analysis, treating 
the liquid as a one-dimensional semi-infinite solid with conduc- 
tion heat transfer. It is believed that evaporation taking place 
at the meniscus at the liquid-vapor-wall interline, neglected in 
all analyses to date, can play a significant role relative to 
that occuring at the bulk liquid interface [10]. 

II. ANALYSIS 

A schematic of the venting system is shown in Figure 1. The 
venting of vapor only is considered here. It is assumed the 
P v (t) is sufficiently greater than P« at all times so the venting 
rate itself can be described reasonably well in terms of choked 
flow, assuming that only single phase flow occurs and that an 
appropriate flow coefficient can be assigned. It is further 
assumed that all the heat transfer processes of significance are 
one-dimensional , whether at the inter facial area Aj_ in Figure 1 
or at the meniscu s L m in Figure 2 . At the present time these are 
the only domains in which heat transfer and phase changes will be 
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considered to take place . The effect of microgravity relative to 
other gravities is presumed to change only the absolute and rela- 
tive magnitudes of the interfacial area Aj_ and meniscus length 
L m . The processes occuring in the vapor space due to the 
pressure decrease with venting will be considered first. 


1 . Vapor Region 

As the pressure decreases, the initially saturated vapor at 
state (T) remaining in the container can be presumed to have 
undergone a reversible adiabatic expansion. This ignores, for 
the moment, any interaction with vapor which might be generated 
at the liquid-vapor interface. Referring to Figure 3, it can be 
noted that if local equilibrium exists during the process the 
state will be at (jT) f a temperature decrease will have taken 
place, and liquid must have formed within the vapor to produce 
the quality required by the state. Since the liquid interface 
and any vapor generated at the interface would be at the instan- 
taneous saturation temperature corresponding to the system 
prssure, with equilibrium no temperature gradients would exist 
within the vapor space, nor need they be considered; the temp- 
erature and pressure within the vapor space are coupled. 

However, the presence of quality in the vapor space introduces 
two aspects in the problem not considered heretofore: The 

quality itself must now be included within the state description 
of the vapor , and homogeneous bulk nucleation must take place 
within the vapor space to produce the uniform quality assumed . 
Homogeneous bulk nucleation requires that the vapor itself be 
subcooled to some degree [11]. The container walls could provide 
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nucleation sites, requiring an amount of subcooling so small as 
to be negligible. One limit then is to neglect the conditions 

necessary for nucleation, assuming that thermodynamic equilibrium 
exists in the vapor space at all times. 

If no nucleation sites were available, either in the vapor 
bulk or on the walls, for an expansion to pressure P2 the vapor 
will follow the isentropic process to the metastable state (2M)in 
Figure 3. This is at a lower temperature than the equilibrium 
state ( 2 ^) because no latent heat associated with phase change was 
given up. For this case, the liquid-vapor interface would be at 
temperature T2, the T sa t corresponding to P2, while the bulk 
vapor would be at T2 m» A temperature gradient thus would exist 
within the vapor space, which must then be taken into considera- 
tion because of the associated heat transfer. The other limit to 
that cited above, which still takes the vapor to be at a uni- 
form state, is to consider that the vapor itself exists in a 
metastable state, with no heat transfer interaction with the 
interface, but that instantaneous perfect mixing takes place with 
the vapor produced at the liquid-vapor interface. The resulting 
mixture would itself be at a metastable state. 

To compare the effects of treating the expansion of a satur- 
ated vapor with equilibrium into the quality region with that 
assuming that no condensation occurs, considering that the vapor 
remains a pure vapor in metastable equilibrium (termed a pseudo- 
vapor here), sample computations will be presented below for the 
isentropic expansion of the contents of a tank initially contain- 
ing only vapor . Certain assumptions will be made for convenience 
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here , to be eliminated when returning to the original venting 
problem. 

Assumptions : 


1 

ii 


111 


The initial contents of the tank is saturated vapor only. 
The contents remaining in the tank at any point will be 
considered to have a uniform state and to have expanded 
reversibly and adiabatically . For computation purposes 
here, the tank will be vented until one-half of the 
original mass remains. 

The expansion of the tank's contents will be treated as a 
non-rate process. With a uniform state assumed within the 
tank, the state of the contents is independent of the rate 
of venting. 

In Figure 3, the initial saturated vapor state is designated 
as (^)f and the final specific volume will be V2 = 2vi with one- 
half of the original mass vented. State ( 2 ^ on the T-S diagram 
is the isentropic end state in the equilibrium quality region, 
while state (2 p) is the isentropic end state with pseudo-vapor 
behavior, with no condensation. Note that V2 = V2p, and the 
problem is now to compute and compare the corresponding end 
states P2, T2 and P2P, T 2P* In general, state (2 p) will be dif- 
ferent than state (2 m) , the difference depending on the degree of 
expansion and on the particular fluid properties. Computations 
below will demonstrate this. 

In both of the above processes, the state of the fluid 
vented is assumed to have the same property as the bulk property 
of the fluid in the vessel . For the isentropic expansion in the 
quality region, the generation of quality requires the formation 


- 8 - 



of liquid drop nuclei , in the practical sense . These nuclei will 
form either as a bulk homogeneous nucleation process or a hetero- 
geneous nucleation process taking place on the walls of the con- 
tainer . For the case where bulk liquid is present in the con- 
tainer, such condensation could also occur on the liquid-vapor 
interface. Since bulk homogeneous nucleation requires a sig- 
nificantly greater degree of supersaturation in the vapor 
relative to heterogeneous nucleation on the walls, it can be 
expected that the latter would be much more likely to take place. 
The implication of this to the venting process is that with the 
droplet nucleation and condensation now taking place on the 
walls, the bulk fluid vented would have the state of a saturated 
vapor only, while the state of the contents of the vessel taken 
as a whole would have quality. The analysis to determine the end 
state for a given mass vented would be different from the two 
cases described above, and would give an end state also different 
from these. 

A brief description of the analysis for each of the three 
cases mentioned above will be given below, to be followed by the 
results of sample calculations for R-ll as the working fluid. 

Case A - Isentropic Expansion with Equilibrium Quality 

Referring to Figure 3, the following conditions apply: 

Initial Conditions: T x , P x (P sat ) (1) 

vi = v gl 

Since m 2 = l/2m x , it follows that 

v 2 = 2v x = vf 2 + X 2 Vf g 2 (2) 

= S f 2 + X 2 S f g 2 (3) 


s 2 = S 1 
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Using tables of thermodynamic properties , T2 ( P2 ) and X2 can be 
determined from Equations ( 2 ) and ( 3 ). 


Case B - Isentropic Expansion as a Pseudo-Vapor 

Referring to Figure 3: 

Initial conditions: Same as Equation (1) above. 

Since it is taken that m 2 = l/2mi, 

v 2 = 2v 1 = v 2 s = v 2 p (4) 

The end state 2P is the extension of the constant volume 
line V 2 into the pseudo-vapor domain such that: 

S 2P = Si (5) 

The entropy change along a constant volume line for low pressure 
vapor with no phase change is given by: 

dS = Cv ^ T ^ dT (6) 

T 

Since state (2 eT) is known, T 2 P can be determined from integra- 


tion of Equation 

(6) as: 

T 2S 

/» 


S 2S “ S 2P = 

S 2S ~ £l 

= / C?<T>dT 

(7) 



J T 




t 2P 



P 2 P is then computed from V 2 p and T 2 P using the appropriate 
equation of state. 


Case C - Isentropic Expansion with Equilibrium Quality and 
Venting of Saturated Vapor Only 
Since the state of the fluid being vented in Cases A and B 

above is the same as the bulk fluid state , the final result ob- 
tained is the same regardless of whether a system or control 

volume analysis is used . For the case being considered here 
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where the fluid being vented is a saturated vapor while the bulk 
fluid remaining in the tank has quality, a control volume analy- 
sis must be used . 


Initial conditions : Same as Equation (1) above. 

As in Cases A and B, it will still be taken that m 2 = l/2m]_. 
The final state will be designated by the subscript "2V", while 
the variable intensive properties in the container will have the 
subscript "a", and the state of the fluid vented will have the 
subscript "e". 

The control volume form of the first law of thermodynamics 
for the constant volume insulated container is given by: 

d(m 0 u a ) = -8m e xh e (8) 

Expanding the left side, and noting from continuity that 6m e = 
-dm a , and that h e = hg a . Equation (8) can be written as 


du. 


hga - u < 


dm. 


m f 


(9) 


Integrating from the initial state "1" to the final state "2V" 
with m 0 2 = l/2m al , 

2V 

m 


/ 


^ = In !ii£2 = lnl/2 = -0.69315 

(hga u a) ^al 


( 10 ) 


Note that u a can be expressed as, 

U a = hf - PVf + X(hfg-PVfg) (11) 

and ( hg a -u a ) as 

hga " u a = h f g(l-X) + Pv f + XPvfg (12) 

The quality in equation ( 12 ) is an unknown variable at the mo- 
ment . Since the contents undergo a reversible adiabatic process , 
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the second law of thermodynamics for the control volume can be 
written, similar to Equation (8), as : 


d(m 0 S a ) = 5m e x S e 

Using continuity 6m e = -dm a , and noting that S e = Sgg, 
(13) can be written as: 

dS g = dm g 

Sgcj — Sg mg 

which is similar in form to Equation (9). Note that S 
expressed as 


(13) 

Equation 

(14) 

c can be 


Sg = Sf + XSfg (15) 

and Sgg— Sg as 

Sga - Sg = Sfg(l-X) (16) 

The integrated form of Equation (9) equals the integrated 
form of Equation (14), and the quality X during the expansion 
must be determined to satisfy this. 


Solutions 

Solutions for the venting of an initially saturated vapor 
from an insulated tank were carried out for Freon-11 for Cases A, 
B, and C above, using the tabulated thermodynamic properties, the 
equation of state, and the constant volume specific heat given in 
Reference [12]. The results are given in Table 1 for two dif- 
ferent initial pressures. 

For Case C, Equation (10) could be integrated numerically 
using either the tabulated property values or the equation of 
state . Case C-l in Table 1 was obtained by numerical integration 

of the tabulated properties in [ 12 ] , while Case C-2 is the result 
obtained with the equation of state . In both cases it was neces- 
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sary to determine the unknown quality at each step by trial-and- 
error in order to satisfy Equation (9). 

It may be noted that, although the final pressure for the 
expansion as a pseudo-vapor (Case B) is just slightly less than 
that for the expansion in the isentropic quality region (Case A), 
the final temperature is lower. This is to be anticipated, since 
the internal energy decrease for the expansion with no phase 
change must take place more at the expense of the vapor internal 
energy, with the attendant decrease in temperature. 

For the Case C, where only saturated vapor is vented, the 
final pressure and associated saturation temperature are lower. 
The final quality is slightly lower than for Case A, reflecting 
the fact that any condensed vapor is retained within the con- 
tainer. 

For the modeling of the original venting problem, where the 
initial contents of the insulated tank consists of a mixture of 
saturated liquid and vapor, it appears that Case C above will 
represent most realistically the process actually taking place. 
The isentropic expansion of a saturated vapor results in conden- 
sation. Condensation nucleation sites must be present for this 
to occur, and it is highly unlikely that bulk phase homogeneous 
nucleation will take place with the presence of metal walls un- 
less the tank is very large. However, nucleation and condensa- 
tion will occur preferentially on the walls and on the existing 
liquid-vapor interface, maintaining the bulk as saturated vapor, 
from which the fluid vented is drawn . 

The alternative to the equilibrium treatment of the vapor 
space described here is a non-equilibrium approach, in which 
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spacewise temperature non-uniformities must be incorporated into 
the analysis . Before such complexities are introduced, however , 
it should be demonstrated that such an approach is necessary . 


It is the rate model for the tank venting process that is of 
interest, and before introducing the interaction with the liquid 
the time element will be incorporated in Case C above, when the 
tank contains only saturated vapor, initially, and where only the 
saturated vapor is being vented, the liquid formed due to the 
condensation occuring as a result of the vapor expansion being 
retained within the tank. 


The continuity equation for the contents of the control 


volume consisting of the constant volume tank is: 


dm a = _ Sm e 

dt dt 



(17) 


The rate form of the energy equation corresponding to Equation 
(8) is: 

j * 

§£ (m a u a ) = " m e h e (18) 


Expanding the left side of Equation (18), substituting Equation 
(17) and rearranging, noting that only saturated vapor is 
leaving : 

du_ 

m o — ( hga“ Ug ) ( — m e ) (19) 


The mass flow rate being vented, m e , is modeled using a 
classical choked flow analysis [13]. In the application to the 
vacuum of space , a choked flow assumption can be expected to be 

valid , and the exit mass flow rate will be a function of the 
upstream vapor properties only : 
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( 20 ) 


m = ?V C D A i K D 
e TrtT]T72 


where 


K D = 


A/2 


Y + 1 


K, *V 

1 v /2 


Y + 1 
2(Y"1) 


and 


Kh 


C D A i K D 

Rl/2 


( 21 ) 


( 22 ) 


Solving Equation (19) for dt and substituting for m e from Equa- 
tion ( 20) : 


dt = 


-m. 


m 1/2 ^ 

T V du a 


(23) 


(hgCF Up) K^Py 

Equation (23) can be integrated numerically to find the time re- 
quired to change from the initial state to any other state at a 
lower pressure. T v and P v are related by the vapor pressure 
curve or the Clausius-Clapyron equation, while the expressions 
involving u a and hg 0 are given by Equations (11) and (12). m a 
comes from: 

= V a /v a (24) 

and 

v a = Vf + Xvfg (25) 

During the expansion process, the quality X must be solved at 
each computational step by trial-and-er ror such that Equation (9) 
is satisfied at each state along the path. Temperature is used 
as the computational variable on the right hand side of Equation 
(23). As will be described later , a number of different sizes 
and series of temperature steps were investigated. 
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2 . Liquid Region 


As the pressure in the vapor domain decreases as a result of 

the venting process , the temperature at the liquid interface, Tj_ 
in Figure 1, will decrease corresponding to the saturation temp- 
erature. If the liquid was intially at a uniform temperature 
this means that a temperature gradient will be established in the 
liquid, resulting in evaporation at the liquid-vapor interface. 

As far as the vapor domain control volume is concerned, this 
represents a vapor mass addition to the vapor space as a satur- 
ated vapor corresponding to the system pressure. 

The rate form of the continuity equation corresponding to 
Equation (17) for the vapor space as a control volume is now: 

dm 

_5 = mi - m e (26) 

• 

where mj_ is the mass flow rate of vapor entering the control 
volume as a result of evaporation at the liquid-vapor interface. 

The rate form of the energy equation corresponding to 
Equation (18) for the vapor space control volume is: 


d 

dt 


(% u o) 


m i^i m e h e 


(27) 


Since the vapor entering the vapor space control volume as a re- 
sult of evaporation is a saturated vapor corresponding to the 
system pressure, it has the same state as the vapor being vented, 
designated as hg a in Equation ( 19 ) . Expanding the left hand side 
of Equation ( 27 ) and rearranging : 

du • • 

mgr — — (hgg - Ug)(mj_ _ ITl0 ) ( 28 ) 
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The rate of vapor generation, m^, is determined from the 
conservation of energy equation applied to the liquid-vapor 

interface . Assuming no heat transfer to the vapor , all energy 
transferred to the interface by conduction in the liquid results 
in vaporization of liquid at the interface, given by: 

qi = m ihfg (29) 

For relatively short periods, where the temperature boundary 
layer is small compared to any radii of curvature present at the 
interface, the liquid may be treated as a semi-infinite planar 
solid. The surface area term, A^, will be that corresponding to 
the shape the interface takes in zero gravity. The one dimen- 
sional form of Fourier's conduction equation for the liquid at 
the interface is: 

91 - -^i(g)| x=0 (30) 

Combining Equations (29) and (30) gives 
• -k-. A i (dt) | 

mi = 1 1 \dx/ 1 x ~° (31) 

hf g 

The problem of determining the interfacial mass transfer is 
reduced to determining the temperature gradient of the liquid at 
the interface, which requires that the transient temperature 
distribution in the liquid near the 1-v interface be determined. 
If the liquid near the 1-v interface can be considered to approx- 
imate a one-dimensional semi-infinite solid in it's thermal 
behavior the analytic solution for a step change in surface 
temperature , in connection with the finite form of Duhammel ' s 
superposition integral , can be used to determine the transient 
temperature distribution in the liquid . The time varying inter- 
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face temperature is taken as the saturation temperature corres- 
ponding to the instantaneous system pressure, which must be 

determined appropriately from the system of governing equations . 

Accordingly, the differential form of the governing equation 

and the initial and boundary conditions for the one-dimensional 

semi-infinite solid, initially at uniform temperature T 0 and with 

a step change in surface temperature to Tj_ are: 

9T = * S 2 T (32) 

3x 2 


T(x,o) = T 0 
T(o,t) = Ti 
T(«,t) = T 0 


(33) 

(34) 

(35) 


The solutions for the temperature distribution and for the temp- 
erature gradient are given by [14]: 


~ T i = erf 


T 0 - Ti 


f*(o,t) 

3x 


= - T o - 


2(at) 1 / 2 

i 


(36) 

(37) 


( irat ) l/ 2 

The interface temperature, being the saturation temperature cor- 
responding to the ullage pressure, will be time varying in the 
present case since the pressure will change as the tank is 
vented. This time varying boundary condition Ti(t) is incor- 
porated into the solution using Duhammel's superposition integral 
[14] in the form: 

t 

d6i(s) 


0 ( x , t ) = 6i(o) • \p ( x , t ) + / \])(x, t-s) 


/ 


ds 


ds 


38) 


Here , 


0 ( x , t ) = T( x, t ) - T 0 

0i(t) = Ti(t) - T 0 


( 39 ) 
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We let 


+ <X,t) E |%tl 

0 i ( fc ) 


( 40 ) 


t ( x , t ) is the unsteady temperature resulting from a stepwise unit 

increase in surface temperature, relative to a uniform initial 
temperature. If the increase is kept at zero until a certain 
time t-s, and at that instant raised to unity and maintained 
constant, the new temperature <j>(x,t) may be expressed in terms of 
( x , t ) as 


4>(x,t) = < 


0 , t < s 
^(x,t-s) , t > s 


(41) 


The solution for 4>(x,t) is given by Equation (36), transformed to 
the form of Equation (40) as: 


\Jj(x,t) = § ( x ' t ) = erf 

0i 


(42) 


2 (at) 1/2 

Solution of the system of equations for the venting problem would 
be performed in discrete time steps, and the discrete form of 
Equation (38) is given by: 

n 

6 ( x , t ) = 0i(o) • i|>(x,t) + I A0i m • \p(x,t-s m ) (43) 

m=l 

where 

A®i-m = 9i( s m) " ®i( s m-l) (44) 

Here, n is the total number of time steps into which the process 
has been divided, m is a running time index, Km<n, and A0i m is 
the incremental change in surface temperature, related to the 
system vapor pressure. 

Since the temperature gradient at the interface is needed to 

compute the interf acial mass transfer rate in Equation (31), this 
can be obtained by differentiating Duhammel ' s superposition inte- 
gral Equation ( 38 ) as : 
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30(0, t) 
3x 


(45) 


t 

= (o)l®(o,t) + / M(o,t-s)£^Ai_^lds 

1 3x J 3x ds 

o 


The discrete form of Equation (45) is given by: 

= 6 i |0 )||(0't) + Ji 4S im ’ |±(o,t-s m > («> 

For a semi-infinite solid with a step increase disturbance: 




1 

[ira( t-s m ) l 1 / 2 


(47) 


Substituting Equation (47) into Equation (46), and noting that 
although for a step initial disturbance that 3\p/3x(o,o) = «, 

©i(o) can be taken as small as desired, so that Equation (46) 
becomes : 

30(o, t) = E (48) 

m=l [ira( t-s m ) J 1 / 2 


Equation (48) is the form used to compute the interfacial 
mass transfer rate below. However, another procedure was devel- 
oped initially to approximate the temperature gradient at the 
interface and will be described here for the sake of complete- 
ness . 

The procedure is to compute the instantaneous temperatures 
at a finite number of points in the liquid near the interface, 
using Equation (43), and fit these points to a third order 
polynomial using a least squares fit. The polynomial is of the 
form: 

T = A + Bx + Cx 2 + Dx 3 (49) 

The temperature gradient of the liquid at the 1- v interface, 

x = 0 , is then : 
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B 


( 50 ) 


The number and spacing of the nodes at which the temperatures of 
the liquid are to be calculated, and with which the coefficients 

A, B, C, and D in Equation (49) will be determined, must next be 
specified. Six nodes were taken arbitrarily as being sufficient 
to obtain the four coefficients in Equation (49). Intuitively, 
nodes nearest to the 1-v interface will give the most accurate 
value of the liquid temperature gradient at the 1-v interface. 
The method used was to estimate a temperature penetration depth, 
<$, taken here to be the depth at which the dimensionless temp- 
erature change computed by equation (36) is 95% of the distur- 
bance, or 


0.95 


and 


erf 


5 

2 ( at ) l/ 2 


(51) 


6 = 1.39 x 2(at) 1 / 2 (52) 

The actual penetration depth would be somewhat less than this 
value, since the actual system does not undergo a single step 
change in surface temperature, but rather a transient change in 
surface temperature. The six equally spaced nodes are taken to 
be within the 10% of this penetration depth nearest the 1-v 
interface . 

Now that the temperature of the liquid at each of the six 
nodes near the 1-v interface is known, the constants A, B, C, and 
D of Equation ( 49 ) may be determined . A least squares algorithm 

was used [ 15 ] which determines the polynomial coefficients which 
minimize the error between the data points and the polynomial . 
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A test program was devised to evaluate the effect of the 
fraction of penetration depth used when fitting a polynomial by 
computing the accuracy of the polynomial in predicting the 
temperature gradient at the 1-v interface. The temperature and 
temperature gradient obtained with the above procedure are com- 
pared with the analytical values for a single step change in 
surface temperature, being the most severe test possible. 

Figure 4 gives the relative errors in the surface tempera- 
ture gradient while Figure 5 gives the relative errors in the 
surface temperature itself, as a function of the fraction of the 
penetration depth within which the six equally spaced nodes are 
located. These computations were also carried out with the 
first, second, and third order curves of Equation (49). With the 
nodes contained within a region of 10% of the penetration depth 
from the surface and using a third order polynomial, the error in 
the temperature gradient at the interface is less than 0.5%. 

3 . Combined Liquid-Vapor Region 

The basic equation used to solve for the transient states 
within the tank venting vapor only and containing liquid and 
vapor initially at a uniform temperature and pressure is given by 
Equation (28). The states of the vapor within the tank and the 
liquid arising from condensation from the vapor state are always 
at a uniform temperature, whereas the temperature distribution 
within the liquid will vary in accordance with the solution for a 
semi-infinite solid, as described in the previous section. The 
system pressure is taken as uniform but time varying. 
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Solving Equation ( 28 ) for the time step dt , and substituting 


for m e from Equation ( 20 ) 


dt = 


( v- u °> r;.-K r 


( 53 ) 


m-i is determined from Equations (31) and (48), expressions in- 
volving u a and hg a are given by Equations (11) and (12), m a is 
given by Equations (24) and (25), and P v can be expressed in 
terms of the T v with the Clausius-Clapyron equation or the vapor 
pressure curve. Temperature of the vapor (and hence the liquid- 
vapor interface) is used as the computational variable on the 
right hand side of Equation (53), with time being the dependant 
variable to be integrated. The quality X, resulting from con- 
densation in the vapor space during the expansion process, must, 
be solved at each computational step by tr ial-and-error such that 
Equation (9) is satisfied at each state along the path. It is 
this element of the solution process that constitutes a major 
portion of the computational time. 


III. RESULTS 

The flow sheet is given in Appendix A, while the symbol and 
program listings are given in Appendices B and C, respectively. 

For purposes of comparisons of the computational results 
using different temperature steps in the expansion process , it is 
necessary to use specific input parameters at this stage . The 


sure at the end of the three second vent period and that pre- 
dicted by the analysis developed here . 

Table 2 lists the experimental parameters from [9], along 
with the final pressure measured after venting for three seconds. 

The acrylic plastic cylindrical container was 6 cm inside dia- 
meter by 10 cm inside length, and had a zero contact angle with 
R-ll. Under the microgravity conditions during free fall suffi- 
cient time was allowed for the interface to achieve its equili- 
brium hemispherical shape before venting was initiated. This 
serves to define the surface area from which evaporation can take 
place. Also listed in Table 2 are the final pressures computed 
by the analysis in [9] and by the analysis presented here. These 
differ in how the condensation from the vapor space is treated, 
and in how the mass transfer at the liquid-vapor interface is 
computed. 

It is noted in Table 2 that although minor variations exist 
in the degree of initial filling and in the initial states, the 
major differences between the test runs is that the venting ori- 
fice diameter increases progressively, which corresponds to an 
increase in the venting rate expressed in the fourth column as 
computed initial ullage volumes per second. The tabulated 
measured and analytic relative pressure drops are the pressure 
decreases taking place in three seconds divided by the initial 
system pressure . 

The parameters of system pressure , the corresponding satur- 
ation temperature , vapor space quality, and mass rates being 
vented and evaporated at the liquid-vapor interface are plotted 

in Figures 6-10 , corresponding to test numbers 1-5 in Table 2 , as 
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a function of time up to a three second maximum period, using 

constant step changes in vapor state temperature of 0.5°F for the 

integration of Equation ( 53 ) . 

Results obtained with different degrees and sequences of 

step changes in temperature during the expansion process are 
presented in Tables 3-7, corresponding to test numbers 1-5 in 
Table 2. All computational results presented here were performed 
with the Harris-800 computer, which has 6 megabytes of virtual 
address space, handles 10 6 instructions/second with 24 bits/word 
and a floating point processor. The objective of varying the 
computational step changes of temperature was to examine its 
influence on the results. The smallest constant step used. 
Computer Runs B with 0.05°F, give the largest final pressure for 
all cases. This is not necessarily the most accurate result, 
because of the accumulation of errors associated with the long 
computational period required for this case. The largest compu- 
tational step used, computer Runs A with 0.5°F, appear to give 
satisfactory results and were those used to make the plots of 
Figures 6-10. 


IV. DISCUSSION 


It is noted in Figures 6-10 that as the initial venting rate 
increases , as a result of the larger orifice sizes , that the 
pressure decreases more rapidly, that the evaporation rate in- 
creases correspondingly because of the associated decrease in the 
liquid-vapor interface temperature , and that the evaporation and 
venting rates approach each other . In Run Number 5 ( Figure 10 ) 
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the evaporation rate reaches a maximum and then begins to de- 
crease, as a result of the lower rate of pressure decrease , which 

in turn affects the temperature gradient in the liquid at the 
liquid-vapor interface. 

From Table 2 it is noted that the final pressure computed 
with the procedure developed here gives somewhat higher pressures 
than in the analytic results of [9]. This is believed to repre- 
sent a more accurate result, because of the more accurate compu- 
tation of the transient temperature gradient at the liquid-vapor 
interface. Both computed final pressures are significantly lower 
than the experimental measured ones. This is attributed to ne- 
glecting, in the analysis, the evaporation taking place at the 
meniscus shown in Figure 2, which can be considerable. Addi- 
tional evaporation would tend to reduce the pressure decrease 
rate. 


V. FUTURE WORK 

1. The results presented here should be placed in a generalized 
form so as to be more generally useful. 

2. The analysis should be extended to include evaporation taking 
place at the solid-vapor-liquid contact line. 

3. Investigations should be initiated on the parameters which 
govern the limits on the pressure decrease rate so as to 
prevent nucleation. For example, it is noted in Figure 10 
that the saturation temperature changes by about 73°F in 
three seconds . This means that if the container were suf- 
ficiently large and contained a significant amount of liquid, 
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the liquid would become superheated by this 73°F, and most 
certainly would nucleate, with local boiling taking place 
most likely on the walls. 

4. Experiments should be conducted with large size containers, 
to study effects of wall heat capacities, interface areas, 
and to corroborate the influence of meniscus evaporation. 


VI . NOMENCLATURE 

a Thermal Diffusivity 

A Area 

Cd Flow Coefficient 

o 

C v Ideal Gas (low pressure) Specific Heat at Constant 

Volume 

d Differential 

h Enthalpy 

hfg Latent Heat of Vaporization 

k Thermal Conductivity 

K Constant Defined Locally 

Kd see Equation (22) 

m Mass 

P Pressure 

q Heat Transfer Rate 

R Gas Constant 

s Time as Integration Variable - Equation (38) 

S Entropy 

t Time 

T Temperature 
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u Specific Internal Energy 

v Specific Volume 

x Coordinate 

X Quality 

y Ratio of Specific Heats 

6 Boundary Layer Thickness, Differential 

0 Relative Temperature 

<j> Dimensionless Temperature 

Dimensionless Temperature Solution for Unit Temperature 
Disturbance 
Subscripts : 

e Exiting Control Volume 

f Liquid 

g Vapor 

1 Entering Control Volume 

v Vapor 

a Control Volume 

sat Saturated Conditions 
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Fig. 1. 


Schematic of venting system. 
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Fig. 7 


Computed results for Run No. 2 in Table 2 


MASS RATE - lbm/sec 










TIME - SECONDS 


Fig. 9. Computed results for Run No. 4 in Table 2. 
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MASS RATE - lbm/sec 


PRESSURE - PSIA VAPOR SPACE QUALITY 



Fig. 10. Computed results for Run No. 5 in Table 2. 
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Table 1 . Computation of final state for venting Freon-11 to 
one-half the initial mass . 


Initial State 

Case 

Final State 

T(°F) 

P(psia) 

X 

T = 7 4 . 000°F 

A 

38.300 

6.754 

.9854 

P = 14.447 psia 

B 

29.138 

6.715 


X = 1.0000 

C-2 

37.989 

6.705 

.9792 

T = 160 . 000°F 

A 

112.208 

28.938 

.9943 

P = 60.451 psia 

B 

109.202 

28.913 


X = 1.0000 

C-l 

112.028 

28.864 

.9920 


C-2 

112.085 

28.892 

.9928 
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Table 2. Experimental Parameters from Reference 9. 

Vessel: Acrylic Plastic - 2.36 in (6 cm) I.D. x 3.94 in (10 cm) long 

Fluid: R-ll (CCi^F) 

Inside Volume: 17.33 in^ (2.84 x 10 ^ m^) 

Venting Time: 3 seconds 


Test 
No . 

% 

Initial 

Vapor 

Volume 

Nozzle 
Dia . 
In. 

Discharge 
Coeff . 

C D 

Reduced 
Flow Rate 
Ullage/Sec 

Initial 

Pressure 

psia 

Initial 
Temp . 
°F 

Final 

Measured 

Pressure 

Drop 

Relative 

Measured 

Pressure 

Drop 

Final 
Analytic 
Pressure 
(Ref. 5) 
psia 

Relative 

Analytic 

Pressure 

Drop 

Final 

Pressure 

Computed 

Here 

Figs. 6-10 
psia 

i 

68 

0.016 

0.64 

0.035 

13.00 

70.07 

12.50 

0.04 

11.83 

0.09 

11.85 

2 

71 

0.035 

0.69 

0.17 

12.75 

70.79 

10.20 

0.20 

8.17 

0.36 

8.36 

3 

67 

0.042 

0.86 

0.33 

13.20 

68.99 

8.80 

0.33 

5.90 

0.55 

6.27 

4 

68 

0.052 

0.875 

0.51 

14.10 

74.03 

7.80 

0.45 

4.26 

0.70 

4.82 

5 

68 

0.076 

0.77 

. 

1.12 

14.65 

72.05 

6.00 

0.59 

1.90 

0.87 

2.71 



Tabl 


Computer 

Run No. 


Constant DT = 0.5 F 11.853 

See Figure 6 

Constant DT = 0.05 F 11.905 

9 Steps DT = 0.05 F 11.838 

then remaining DT = 1 F 

Geometric Series Init. 11.890 
DT = . 025F Ratio = 1.05 

Same as D, but 11.890 

DT max = 1 • 0 F 

Geometric Series Init. 11.886 
DT = .05 F Ratio = 1.05 

Same as F, but 11.886 

DT max = 1 • 0 F 


al Vali 

aes at 3 £ 

3ec . 

X 

# 

m vent 

m evap 

.998 

. 719E -4 

. 230E -4 

.998 

. 722E“ 4 

. 256E~ 4 

.998 

. 718E -4 

. 214E -4 

.998 

. 721E -4 

. 241E -4 

.998 

. 721E** 4 

. 241E -4 

.998 

. 721E -4 

. 240E -4 

.998 

. 721E -4 

. 240E -4 








ble 4. Computer Runs for Test Number 2 in Reference 9 


Computer 
Run No. 

Temperature 

Steps 

Final Values at 3 Sec . 

P 

X 

m vent 

ra evap 

2k 

Constant DT = 0.5 F 
See Figure 7 

8.361 

.990 

. 266E -3 

. 111E -3 

2B 

Constant DT = 0.05 F 

8.454 

.990 

. 269E -3 

. 116E -3 

2C 

9 Steps DT = 0.05 F 
then remaining DT = 1 F 

8.312 

.990 

. 265E -3 

. 107E” 3 

2D 

Geometric Series Init. 
DT = . 025F Ratio = 1.07 

8.340 

.990 

. 266E -3 

. 106E~ 3 

2E 

Same as D, but 
DT max = 1 • 0 F 

8.345 

.990 

. 266E -3 

.107E“ 3 

2F 

Geometric Series Init. 
DT = .05 F Ratio =1.07 

8.337 

.990 

. 265E -3 i 

. 106E -3 

2G 

Same as F, but 
DT max = 1 • 0 F 

8.342 

.990 

. 266E -3 

. 107E -3 





5 • Computer Runs for Test Number 3 in Reference 


Computer 

Run No. 


3A 

3B 

3C 


Temperature 

Steps 


Constant DT = 0.5 F 
See Figure 8 

Constant DT = 0.06 F 

9 Steps DT = 0.06 F 
then remaining DT = 1 F 


Final Values at 3 Sec, 


6.269 

6.364 

6.213 


X 


.979 

.980 

.979 


m vent 


-3 


365E 


370E" 3 


362E 


-3 


m 


evap 


. 187E~ 3 

. 194E -3 
. 184E -3 


3D 

3E 

3F 

3G 


Geometric Series 
DT = .03 F Ratio 

Same as D, but 
DT max = 1 • 0 F 

Geometric Series 
DT = .06 F Ratio 

Same as F, but 
DT max =1*0 F 


Init . 

= 1.06 


Ini t . 

= 1.06 


6.203 

6.232 

6.203 

6.230 


.979 

.979 

.979 

.979 


. 361E 
. 363E 
. 361E 
. 363E 


-3 

-3 

-3 

-3 


. 178E 
. 183E 
. 178E 
. 183E 


-3 

-3 

-3 

-3 








Table 6. Computer Runs for Test Number 4 in Reference 9 


Computer 
Run No. 

Temperature 

Steps 

Final Values at 3 Sec. 

P 

X 

• 

m vent 

• 

m evap 

4A 

Constant DT = 0.5 F 
See Figure 9 

4.822 

.966 

. 439E -3 

. 258E -3 

4B 

Constant DT = 0.06 F 

4.909 

.966 

. 447E -3 

. 266E -3 

4C 

9 Steps DT = 0.06 F 
then remaining DT = 1 F 

4.771 

.965 

. 435E -3 

. 254E -3 

4D 

Geometric Series Init. 
DT = .03 F Ratio = 1.04 

4.757 

.965 

. 434E -3 

. 249E -3 

4E 

Same as D, but 
DT max = 1 • 0 F 

4.784 

.965 

. 436E -3 

. 254E -3 

4F 

Geometric Series Init. 
DT = .06 F Ratio = 1.04 

4.757 

.965 

. 434E -3 

. 248E -3 

4G 

Same as F, but 
DT max = 1 • 0 F 

4.783 

.965 

. 436E -3 

. 254E -3 
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Table 7 . Computer Runs for Test Number 5 in Reference 9 


Computer 
Run No. 

Temperature 

Steps 

Final Values at 3 Sec . 

P 

■ 


« 

m evap 

5A 

Constant DT = 0.5 F 
See Figure 10 

2.711 

.927 

. 476E -3 

.356E" 3 

5B 

Constant DT = 0.05 F 

2.771 

.928 

. 486E -3 

. 364E -3 

5C 

9 Steps DT = 0.05 F 
then remaining DT = 1 F 

2.677 

.926 

.470E" 3 

. 351E -3 

5D 

Geometric Series Init. 

DT = . 025F Ratio = 1.028 

2.656 

.925 

. 467E -3 

. 346E -3 

5E 

Same as D, but 
DT max = 1 • 0 F 

2.682 

.926 

. 471E -3 

. 352E -3 

5F 

Geometric Series Init. 

DT = .05 F Ratio = 1.028 

2.657 

.925 

. 467E -3 

. 346E -3 


5G 

Same as 

F, but 

2.682 

.926 

. 471E -3 


DT max = 

1.0 F 





. 352E 


-3 
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Appendix A. Algorithm Flow Chart 
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Appendix B. Symbol Listing 


FORTRAN 



Symbol 

Description 

Units 


AIN 

Total mass of evaporated vapor 

Lbm 

AOUT 

Total mass of vented vapor 

Lbm 

AMASS 

Initial mass of vapor domain 

Lbm 

AS 

Interface surface area 

in 2 

AT 

Nozzle corss sectional area 

In 2 

CD 

Discharge coefficient 

- 

COND 

Thermal conductivity of liquid @T 

Btu/hr-ft-°F 

CV 

Constant volume spceciflc heat of vapor @T 

Btu/Lbm °F 

CP 

Constant pressure specific heat of vapor @T 

Btu/Lbm °F 

DOTM 

Mass flow rate across interface 

Lbm/Sec 

DT,DTEMP 

Temperature step 

°F 

DTI ME 

Time step 

Seconds 

DTR 

Temperature ratio 

- 

HF 

Enthalpy of saturated liquid @T 

Btu/Lbm 

HG 

Enthalpy of saturated vapor @T 

Btu/Lbm 

OUT 

Mass flow rate of vented vapor 

Lbm/Sec 

PSAT 

Saturation pressure @T 

Psia 

SF 

Entropy of saturated liquid eT 

Btu/Lbm °F 

SG 

Entropy of saturated vapor @T 

Btu/Lbm °F 

T 

Temperature 

°R 

TF 

Final ullage temperature at next step 

°R 

Tl 

Initial ullage temperature 

°R 

TIME 

Total elapsed time 

Seconds 

TIME2 

Elapsed time for each step 

Seconds 

TO 

Initial ullage temperature 

°F 

UF 

Internal energy of saturated liquid @T 

Btu/Lbm 

UG 

Internal energy of saturated vapor @T 

Btu/Lbm 

VF 

Specific volume of saturated liquid @T 

Ft 3 /Lbm 

VG 

Specific volume of saturated vapor @T 

Ft 3 /Lbm 

XF 

Final quality of vapor at next step 

- 

XI 

Initial quality of vapor 

- 
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Appendix C. Program Listing 


mil 

•W-U 

Pi§§p 

'fej 

.. .XtSC 

C 

c 

c 

n 

c 


B I HENS Li IN DT ! .*£ ( 1 OOO ) , DTEMP < 1 OOO } , COB (5) , ATT (5 ) , AM ABBS ( 
u TOO ; 5 / . DTT ( 5 > , DTRR ( 5 > 
v 1. ai'-H AKBiE uOErrj LIEN'’" 

;>AT 000/ ,o4,.6V, .86, .875, . 77/ ' . 

ARE A iT" Uliil i 'i : 2*2 £. A i THEGAT If'] OMUAn’c I N 

DATA AT i / 2 . C0666E-4 ,9.621 IE-4 ,1.3 9376E--3 ,2. 12ii4iE-3, 

ft i :T 5 4 5 1 1 i - 3 /' 

iNXL IAL VAPOR iiil IN LBM ’ 

DATA AHA3S8/2, 225241 iE-3,2. 2765154E-3',2. 2225I81E-3, 

&2. 4005 C ?19E~3,2.4879326E~3/ 

INITIAL ULLAGE -TEMPERATURE IN DEG F 

DATA TOO/ £tS ... 67 95 , 67 . 730 , 69 . 45 1 , 72 . 760 , 74 . 707/ 

INITIAL TEMPERATURE STEP IN DEG F 
DATA DTT/. 05,, 05,. 06,. 06, .05/ 

TEMPERATURE RATIO 

DA VP. DTRR/ 1 . 05 , 1 . 07 , 1 . 06 , 1 . 04 , 1 . 028/ 

LI QUID- VAPOR INTERFACE AREA IN SQUARE IN 
AS—8 . 6303608 

RUN 5 TIMES BASED ON 5 EXPERIMENTAL. DATA (INPUT) 

DO 1 L“1 ,5 

PRINT RUN NUMBER , INITIAL TEMP, TEMP RATIO 
T RITE (20 n 500) 

m 1 -€ ( 20 , 600 > L , DTT (L) , DTRR CL ): 

'SET iTlTfAL QUALITY OF VAPOR DOMAIN 

XI-1. '. ' 

SET DISCHARGE COEFFICIENT 

CD" ODD (L) eS ■ ' ; 

ST r AREA OF VENT NOZZLE AT THROAT 

SLV INITIAL VAPOR MASS 
AHASS-AMA8SS <L) 

SET INITIAL ULLAGE TEMPERATURE 
TO :0 ( L 

SET INITIAL 1 EMPLRATURE: 'STEP 
DT— DTT (L) 

SET TEMPERATURE RATIO 

wmtmmm 

PRINT HEADING 

WRITE (20, 100) 



INITIAL TEMP FROM FAHRENHEIT TO RANKINE 

V. T Z'/’i * ’’ L-TS 


TOTAL ELAPSED TIME 

TIME-0. ' 

INITIAL MASS FLOW RATE ACROSS LIQUID-VAPOR INTERFACE 
DOTH 1-0. 

TOTAL MASS OF VAPOR ACROSS LIQUID- VAPOR INTERFACE 
AIN— 0. ' ' '" : - 

TOTAL MASS OF VENTED VAPOR 
AOUT-O. 


51 



St:: 

*; 

.&] 

£• 

C 

r'' 

C- 

c 

c 


c 

c 


SPRINT ^ 

TIfNE=TOTAL ELAPSED TIME 11 SECOND 

11 11 1 iilMiiiii 1 1 H Hill 1 1 jpill i— « 

MHnfellli HHHK PREBSIjEE IN Ka 

f|§f§pll nv OF VAPOR DOMAIN ipOgpf 

■' AA'oar^siIsiiipbo 

DOTi'li— HAS3 FLOW PATE ACROSS LIQUID --VAPOR INTERFACE IN 
007= MASS FLOW RATE OF VENTED VAPOR , N LuM/SEC 

A [M- i OTAL MASS OF VENTED VAPOR IN LBM 

AOUT=TOTAL MASS OF VAPOR ACROSS LIQUID-VAPOR INTERFACE 

WR I TE ( 20 , 200 ) T I ME , TO , PEAT <TI> , X I , AMASS . DOT M i , OUT ( T : 

T 5 A I N , AOUT 

SET ALL TIME STEPS AND TEMPERATURE STEPS ZERO 

DO 6 J*i , iOOO ; : ' • 

DTI ME < j ) — 0 . 

DTE iV T- LT ) =0 . 

6 CONTINUE ' ■ ■ 1 ! ' 


2 CALL STAT ( T I - X I „ DT , TF , XF , A ) ■ . 

j HMSl® .. . 

DT I ME ( K > “VENT < T 2 , X I , TF , XF , CD , AT , AMASS , A* AIN ASS ) 

T I ME 1 * T I ME +DT I INE ( K ) 

DTEM;- ! (K> =DT 

4 DO 5 1=2 , K 

T I ME. 1»TIMEJ - DTI ME (I - 1 ) 

’ DOTM2*DOTH24-eVAP <TF , DTEMP < I > , AS **/&!£* K . : - i ' ,-/;l 

5 CONTINUE 

T I M£2*fiVENT ( T J , X I 4 TF , XF ,CD , AT v AM • ;8S , A* AMASS , DQTM t t 

7 0L---TIME2*,! 0000 i 

: F < : )BS i T I ME j. - ", I i IE 2 ) . LT « TOL > T HE N 
SC ! U > 

ELSE 

DSliii® ■ . A ; ■- : 

T 3 ME 1 “ T I MET’- T | ME 2 
80 TO 4 

END IF I | : : I II III! I II II I 

3 DTI ME CIO r-TIME 2 


Ui 1ME2 

AIN—AIN+ <D0TMH-D0TM2> *. 5*DTIME (K) , ’ I . ^ : ; T J 

iMIntt (OUT a I SB , AT > +QUT < TF .CD .AT) )*..B*DT 2 ME £ id 



TI=TF 
X I = XF 

D0TM1=D0TM2 7 Ib 

K~K‘*l : B ; TAVa 

AMASS*A*AMASS 
■DTLDT *DTR ; 

IF(i;- IN 1, ) THEN 


i 

IOC 

!~!UO 

-2'i z J'k? 

: ,, CM:v| 

fell; 

!i!iBllit 


t 


' ' " 1 ‘ R 020 200 Of THE IP---459. 67 , PEAT ; IF) .XF .AMASS, DOTH! , 

■CUT < 7 F. CD, AT, .AIL. A,, Of . ' 1 ' 

!F (TIME. ST. 3.0) GO T L.t I 

continue 

.-•JKMai.;// , 4X, 'TilH-. ’ ,4X, 'TSAT ' ,SX , ' PSAT ' ,5X , 'QTY ' ,6* . 'VAP MASS' 
I RATE ' , 4 VENT RATE' ,3X, 'MASS EVAP ' ,SX , ' MASS VENT' ,/) 

FDkMAT <2X ,F8. g:,;tX ,FS. 4 , IX ,P9. 5 , 1 X , r 7 - 5 , S ( 2X , El 2 Jo ) ; 

FORMAT 4X , ' <SEC) ' , 3X.., : .' <DE8 F) ’ 1 2X , • (PSTA) ' , 12X , ' (LSM> V .9X, 
^LBM/SEC) ' ? 5>; , ' v-LBM/SEC) ' . / } 

FORMAT (./// / „ 4X , ' RUN #',7X,'INIT TSTEP ' , 6X , ' TRATl O' ) 

FORMAT C6X s ll , 10X,F6.4, JOX,F3.3) 

0jfj» ■ : 

EM) ''''' " : ' WWIM 


C : 
C 


; M FTC T i ON COhiD CALCULATES THERMAL CONDUCTIVITY GF LIQUID 
MMCN-li AS A FUNCTION OF TEMPERATURE 
.MUDLCn V1TY 13 IN EMU/ (HR) (FT: (DEG F) 
rSiMMAATURE IS IN DEG RANK. I NE 

FUNCTION COHD(T) . ; . 

I F C f LT . 4 1 9 67 ■>. 1 HEM 

C. ! .iND*i 0633 ' ] 

EL.rE 1 F ( r . LT .. ASF . u7 > THEN 

Ci iNU 35 . 0633+ ( . OUST - 0633 ) * < 1-4 i 9., 67 ) /40. 

Hi. bir 1 Ft 1 , LT 499 . 6 7 ) THEN 

Ctmti-. 0587+ ( . 05 4 C 058/ ) * (T-459. 67) /40. 

EL SEIF <T . LT. S79. 67) THEN 

i.* i. i'V! /— . u 49 — . oL*4o ) 4 i f —479 .67 . 

c 044 9 : 

HKD IF 
RETURN 

s ENl! . ... 


,iu 

g|§ 

c 

c 

|§S 

m 

c 

c 

c 


FUNGI ION EVAP CALCULATES MASS FLOW RATE ACROSS LI QUID- VAPOR 
INTERFACE DETERMINED FROM CONSERVATION OF ENERGY EQUATION 
APPLIED TO LI QUID- VAPOR INTERFACE 
EVAP is IN LBM/SEC 

T IS TO CALCULATE PROPERTIES IN DEG RANK 1 NE 
DEI..T IS TEMPERATURE DROP AT LIQUID— VAPOR INTERFACE 
AS IS LIQUID-VAPOR INTERFACE AREA IN IN2 


FUNCTI ON EVAP IT , DELT , AS, TIME) 
EVAF*-SQRT ( CONI) < T > * CP ( T ) / < T I ME 
£■/ VMfHT) -HP < T ) ) / 1 #314 .00127 


) ) *AS*D€LT 


c 

as 


FUNCTION llll CALCULATES TIME WITH EVAFGRA* I ON TAKER I17NEEN 
2 , ■ 



* ( US ( T 2 , X 2 ) -U 0 ( T 1 , X I ) J 
RETURN 

s® 




W ;: 

If* : 

§§.) 

c 


i function i sate 

^^^^pripN.; FOR 

FUNGI I ON RATE ( T , X , CD , AT , AMASS ) 

hMHK$ li? 5 » i 

H^l RETURN 



aT,X>*PSAT (- n* 2. 


c 

kr*f 


FUNCTION OUT .CALCULATES MASS FLOW RATE OF VENTED VAFC- 
OUT IS IN UBH/3EC . •■ ■ ' 


T IS IN DEB RANKING 

**" ImMHHHR*' 


I IN ; .;■■■: 


ISIiliil 


FUNCT 1 DM HI IT i T .. PT> .. A 1 


FUNGI 1 ON OUT i T CD , AT 

1 1 6492 7 * Ifpgp 


• ■ ' " • * / w JL lOV7^./ 

CUT =?'£ AT (T) *CD*AT*AKD(.T> 7 (SORT CT*R> *2. 11507642) ' 1 

RETURN • 


C 

c 

II 

III 




r^^ION A VENT TIME W I THtJUf ** ,wr ‘ 


$ . .ww . . .. w <. . , mmmmm 

BE! WEEN 2 GIVEN CONDITIONS 

FUNCT I ON A VENT (T 1 , X 1 , T 2 , X 2 , CD , AT , AMASS , DM ASS , DOT M i . I DTM2 ) 
A VENT- < ARATE <T1 , XI .CD, AT , AMASS , D0TM1 ) 

r*tv * -r*' vs** a #»>/*•« v SfS >i ir * >v 4i '"fit* ><T- * , * * v 


T* ARATE ( T2 , X2 , CD , AT , BMASS , 
RETURN 


C 

ft 

i 

c 




FUNCTION ARATE CALCULATES l/(dU/dt) INSIDE VAPOR DOMAIN 
WITHOUT EVAPORATION FOR GIVEN CONDITION 
FUNCTION ARATE <T,X, CD, AT, AMASS, DOTH) 

R=. 0781 164927 

ARATEw- AMASS /FUN <T , X ) / <CD*AT*AKD <T)*PSAT (T) / (SORT (T*R> * 
H55764t)-liOTM) » ' illll 111111 



END 


. CUD \ ION AKD CALCULATES VAUJE DF' CONSTANT IN CHOKE FLON EQUATION 

llllllll lllllp ■ ■ ■ 111 

gamma~cp co /cv < t> . ... 

AND--SQRT t GAMMA > .* < 2 , / ( GAMMA* 1 .>)**< < GAMMA* i . >*.r. / « GAMMA- 1 ) > 

IIS? Sllillt x T .to, ;;tt 01:10:^0 ;: : t 


•r * : 
1 

Q. 

Ill 

III 


SULRUUiTNE ST AT CALCULATES QUALITY AND VOLUME CHANGE OF VAPOR 
DUMA IN GIVEN INITIAL TEMPERATURE, QUALITY, AND TEMPERATURE 

drop ..I;.- ’ 

SUBRGUT I NE STAT ( T 1 , X 1 , DT , T2 , X2 , A ) 

VOL i —VOL ( T 1 , X 1 ) 

T2-T1-DT 

i* < & < V ^ , X 1 ) -SF < T2 ) ) / • 3G < T2 ) -SF ( 72 ) ) 

VALUE i=D IF <T1 , XI , T2 , X2) 

DX»DT*. OOOOl 
1 a2-L2-DX 

VALUES— Di I F ( T i , X 1 , T2 , X2 1 

< ABS < VALUE 1 ) „ 81 . ASS ( VALUE2 ) ) THEN ' ' 

- 'TALUc.I -VALUES 
GO TO 1 

, >:2^.<2+nx . ; . : ; . 

LOT IF V ; 

A- VOL I / VOL. ( ! 2 , X 2 ) 

RETURN 

END 


i 

III 

I# 


■G y 

C 

IS 


c 

c 


FUMCT I ON FUN CALCULATES VALUE OF (HG-U) INSIDE VAPOR DOMAIN 

- UN» ( HG ( T ) -Hr ( T ) ) * C 1 . - X ) +PSAT < T ) * VOL ( T , X ) * . 1 85053 
RETURN 


FUNCTION UG CALCULATES INTERNAL ENERGY INSIDE VAPOR DOMAIN 

FUNCTION US n , X > 

UG==HG<T)-FUN(T,X) 

RETURN . 

’ END . .■ ... ■: - ..: v . 


-i» 






VOLUM 


function voi_a,x> 

VOL— X *V 6 < T i+!L -X ) *VF ( T ) 
RETURN 

ill Hi i i 1 H Mi i 


VI ON s CALCULATES entropy 


up VA-CR DOMAIN 


FUNCTION S(T,X> 

S- X *SS < T ) 4 - ( 1 . - X ) *SF < T ) 

RETURN 

END 


IS USED IN SUBROUTINE STAJ TO 


FUNCT I ON D IF CTi , X 1 , T2 , X2 ) 

Z ;L- (US < T2 » X2) —US (Ti , Xi > } «. 5* < 1 . /F 
Z2=AL0S ( VOL CTi , XI ) ) - ALOE CVGL (T2 . X 
DIF— ZJ— Z 2 
RETURN 


XI) +1. /FUN <T2,X2)> 


RALCULAt^ 


: iiATED j PRESSURE; 


;#alog: 


T) +D*T->-E* <F-T> /T*< 


<10. >*P1> 


FUNCTION VF CALCULATES VOLUME OF SATURATED LIQUID FREON-11 

vf is in ctf ft/lbm ’ 

T IS IN DES RANK1NE 

FUNCTION VF(T) 

A=34. 57 
B=57 . 638 1 f 






/ 3 . ) +C*F*#.<(2» /3 


•T I ON . ANEWTQN WILL BE USED IN • FUNCTION V3 TO FIND 
"OF EQUATION OF STATE BY USING NEWTON ' S METHOD 

UNCTION ANEWTQN ( T , V ) 


B~. 0019 


■3. 126759 
•00 1 3 18523 
35. 76999 
. 025341 


A5™—2. 3S893E— 5 
fe 5 — ' 2 . 4 4 8 3 0 3 E - 8 
C5-—1 , 47S379E-4 

|6«-9.4/2t03E4 

« O / 

K k ' 4 1 , 5 ' 

H-AK:.Y r T/TC 


■ - ^ " ' B *** R ■** 1 / V P •* \ P 23 •+■ £ 2 & T •+• L’ 2 •* E X F* ( H 5 ) / V B 2 
T < A3+83*T +C3*E XP (H) ) /VB**3. + (A4+B4*T) / 
| (.ASYBillillllll 111 1 /VB**5. 

i.)F~-.R*T/ VB**2. -2, * <A2+B2*T+C2*EXP (H> ) 
43. * >; A3+83*T+C3*EXP (HI ) /VB**4.-4.* ' A4+ 
$5. * <A5+B5*T+C5*EXF* <H) ) /VB**6. 

IF (V.GT..1) GO TO 1 
F=F+ <A6+B6*T> /EXP (580. *V> 

DF-DF-580. * ( A6+B6*T> /EXP (5BC 
% ANEWTON=F/DF 
RETURN 
END 


TT : : > C. 3 ;T; : 


C i;STrt:F 

k FUNCTION VS CALCULATES VOLUME OF SATURATED VAPOR FREON- 11 
C V© IS IN CU FT/LBM 

0 T IS IN DES RANKINE SilT 



1# n n n n o 







r UNCTION VB < T :* 
rF«7-459.67 

IF (IF. LT.-38. ) THEN Jv 
VG®50. 

|| HE IF (TF.LT . — 14. } THEN 

IfiiWililiffi I g : F . i t 111 : 1 g I m ■ 1 il : i! : is ISliil 

| ELSE IF lllllil 13. ) THEN 
VO® 15. 

ELSE I F (TF . LT .43.) THEN 
VO® 7. 

ELSE I F (TF . LT .92.) THEN 
VO =3 . 

ELSE IF (TF.LT . 136. ) THEN 
VO® 1.5 

ELSEIFCTF.LT. 188. ) THEN 
VO®. 7 

ELSE I F ( TF . LT . 266 . ) THEN 

■ . llllilll 

ELSE IF (TF. LT. 200. ) THEN 
VO®. 1.5 

VO®. 11 

.. '.. ; ;.ELSE' gS :£ |J||| 1 1| 

VO® . 06 
END IF 

T OL® VO* 1 E-4 

1 y i *V0-ANEWT ON < T , VO ) 

1 1 ( ASS < V 1 -VO ) . LT . TOL I SO TO 
VO® VI 
SO TO 1 

2 VB®V1 • 

RETURN 




FUNCTION HG CALCULATES ENTHALPY 0 

® §■ t ' I 

TIE I 1 DEG RANK I ME 


ATURATED VAPOR P¥& 


N— 1 1 




FUNCTION HG < T ) 
A® 42. 14702365 
8®— 4344. 343807 
C®- 12. 84596753 
D®. 00400837250 
E= . 03 1 3605356 
F®862. 07 
R®. 0781 17 
BB®. 0019 
A2®-3. 126759 
82®. 001318523 
C2---—35. 76999 



nnno 


. B7C 1 2 .1 E-~3 

■ || S i 1 |p y," O C* *T; / 

h 4- : .. OU l 6c7/ 7 / 
"•• 4 --;; 80 BGe-. 2 E 7~6 
«. 35 d 93 E ~5 

^ 5 ^ 2 , 44 S 3 ^ 3 L — S'" 


®s* i 


,2 


;d= 


OE POOH 'QyWfY 





' ' 1 l-T : ' 


il!i®!!!lll 

111 


v •■ i; ii 

., 7‘575 ( >4 6.3 

ei«~ 9 , 472 i 03 |%;.'' , mgf™ L| v x ^% ; ft|F~.- '\1 

• ;TC“ 84 B. C 7 -l-lf" W :? R 1 "* iV ■ I ia»ll>l»l 

A«*- 4.5 IPII ■ 111 I : " : ' : 1 ! I 

a: =. 0233 15 
Bl=~336.«070 3 
Cl =2, 798S23E-4 
Dl~— 2. 123734E-7 

41 =~S. 777010 4-1 1 .. |j|g . : . 

' T0=-40».+^59,.67 

PT80=PSAT CTO) *ALGG < 10. > * C ~B / TO / T0+C/70 / ALOG ( 10. ) +D- 
$E.*F/T0/TG*AL$S1G CF-TO) -E/70/AL0B < 10. ) ) 

V2-VSCT0) 

V1--7F < Toj •' •: : 

H2~P78o*T0* < V2-V1 > *. 185053 
U2«H2-PSAf *T0)*V2*. 135053' ' 

tiT^AK/TC*T ■ 1 '•'■■■■ 

8 ! 0-- : AK/ i C* ■ 0 : 111 ’ ,': i ||l|f|T 

V. 6 ~ V2--7H . . ' ' ; 2;2 ■ ■■ ■ ■ ■ ' ; v . , 

v r*vs c ' ) 

v~; 2-70 ;T> ~6B - . - ; . ,|1 1 

832= * C?/y3*. 5 * C 3 / v 2 , ' V F+,2 5 * C5 / V & * * 4 . ) * .if 2 ' : ■ 

iiOFS; ) •> (GT— 1 . / - EXP (8 i'O) * ('3TC— :i . > >*, 185053- 
tfAl*<T-TO>+Bi* Cl. /T—i . /TOJ-.S^CI^cT^T-TO-frTOl- 
SDi/S, *fT**3. -T0**3. ) - . 25*E 1 * ( T**4 . ~TQ**4 . > 

U63~ {-A2+C2*(ST- 1 - > *EXP C GT i ) * i 1 . '/VTB- 1 . /VB) + 

4C- -A3+C3*<GT-i. ) *EXP <GT> ) * . 54 ( ,1 . 7VTB**2. -1 . /VB**2. )- 
-$ A 4 * 1 . / 3 .* ( 1 . / VTB **3 .-1. / VB**3 . > + 

s C-A5+C5-MGT-1. > *EXP (GT) ) *. 25* < 1 . /VTB**4.-1 . /VB**4. ) 

2-1)63*. 185053*02 
Sl!f#te * m AT f T ) *V T* . 1 85053 - 

SrrURN - Siipl 0 11 

END 



T IS XN DEB RANK I NE 


FONCTION HF ( T ) 

B* 8 — 4344. 343807 
O- 12. 84596753 
D“ . 00400837 2507 

# - . 03 1 3605356 



U Cl u u o u 


PT8;:-=PSAT ( T) -paLOG < 10 .♦)*< -B/T/T+C/T/ALQS C IO. ) +D- 
*E*F / T/T *AL L-3 1 0 ( F--T) -E/T/ALOG <10. > ) 

HF S-PT Sfa*1 •* (VS 21 ) -OF (T) ) #. 185053 

HF=*H8 I T ) — H r G 

RETURN 


FUNCTION SS CALCULATES ENTROPY OF SATURATES VAPOR FEE 
SO IS IN BTU/LBH i>ES RANK I NE 
T J.S IN DES RANKING. 

FUNCTION SB (T> 

A=42. 14702865 
B— 4344.343007 
C=- 12. 84596753 
0*». 004008372507 
£=.0313605356 
F*862.07 

•|||S|*p3- : 1 1 :■• ^ : 3|3 

88 “-. 0019 

A2=-3. 126759 0 %%' . ■ 

B2». 001318523 

. .. C2*-35. 7699.9 ^ 771 El : |I8!11111 

A3=— » 02534 1 

. B3vls 875121 F.-5 ' 7 .. 3 .:; V;-V ..3 vVfl V|# ^oVV 

C3~ 1 . 220367 

A4=. 001687277 E??t 7 

B4--—1 . 805062E-6 

■ '■■ A5=~2. 35893E— '5 ' ( ■ 7/ • 7. 

85=2. 448303E-8 


mm 


,478379£V4 


Al«i;. 05750418 
•B6»-9.472103E4 
TC-34B.07 
AK---4 . 5 ; .■ 

A! =.023815 
B1 =-336. 80703 , 

Cl =2. 798823E-4 
D'=-2. 123734E— 7 
1 1=5. 99901 BE- 11 
T0=— 40. +459. 67 

PTSO=PSAT < TO ) * ALOB < 10. ) * (-B/TO/TO+C/TO/ ALPS < 10. ) +B 
$E*F/T0/T0*ALQG10 (F— TO) -E/TO/ALOG < 10. ) ) 

V2-VS (TO) 

V1=VF (TO) 

H2=PTS0*T0* (V2-V1 ) *. 185053 
S2=H2/T0 

lliK/il 

8T%iS#^,>- ■ 

GT0=G*TO 


- 60 - 



UCJU POO ; 81 O U O O'! 


VB---V2-E-B 
... VT=M&(T) 

VfB«V8iV>-BB 

S32=8*UI2/VB+.®*C3/V8/VB4.25*C5/VB«*4. >* 

*(eKP<aT)"EKp(S!ro> >*. 18 SGj®3- 

siv-. I * i ALG6 % T ) ~ Ai_OG i PO > /T/T— 1. / fO/1 0) - U.J.* ( I ”1 0) " • 

Lot / 2 . * <T* r~TO*TO > —El / 3 . * <T** 3 . -T 0 ** 3 . > • 

SSliRi lILDBCVlS -ALOS (VTB) > + <B 2 +C 2 * 6 *EXP < 67 ) >*U. /VTB-i . /VB)+ 
4 CB::^C 3 «-Qs-EXP ( 67 ) > *. 5 *< 1 . /VTB /VTB- 1 . /VB/VB) f 
32 : 4 / 3 . *< i, /VTB** 3 .- 1 . /VB** 3 .) + 

•4M8§+C5*8*i~XP (ST) ) 4. 25* (1 . /VTB**4. -1 . /VB**4. ) 

S8*-S32-863*. 185053+82 2' : 

' RETURN 
END 


FUNCTION SF CALCULATES ENTROPY OF SATURATED LIQUID FREON-1 1 
3F IB IN BTU/LBM BEG RANKIiME 


IN DE8 RANKINfc 


FUNCTION SF CD 
B— 4344.343807 
t>- i 2. 84596753 

'D«. 004l*f 72507 ; ? ■ 

E=. 0313605356 ... ■,. 

F=862 . 07 ... /.■T: v 

PtS6*PSAT | |gj #AL06 UC4 ) * (-B/T/T+C/77 AL08 (ID.) +D- 

SFB=PTS6* CVS(T) -VF <T) )*. iS5053 

3F=S8 (T > — SF6 

RETURN 

END ' f v:'' ■ ... 


FUNCTION CV CALCULATES HEAT CAPACITY OF VAPOR FREON- 11 AT CONSTiNT 
VGLUNE ii 

CV It !M BTU/LBR DEG F 

FUNCTION CU(T> 

BB- . 00 1 9 
C2“— 35. 7.6999 

C3“l . 220367 aiiivaia 

C5=— 1 . 47S379E— 4 
Ak=~4.5 ■ 

A 1=. 023815 
B1 =-336:. 80703 ' 

Cl =-2. 798823E-4 
Dl=-2.Ii§#34£-7 
fli Afiloi llBSt 

TC=848. 07 
GK=AK/TC 


* 61 - 



HSSSllPliliiii# 


3'T=i,K*T , 

VT=V8 (T 5 --BB 

c|| (-82^W“C3,/VT/VT-C^/VT^4. 

I £ >; P ( ST ) * 8K ftpllfi ■* T ; * . 1 35033 if A 1 4- 8 §;si 

$1, 1 4f | *& "! 4* E 1 •¥: : , 

RETURN 

END ... : 


'T :i ii 






FUNCTION CP CALCULATES HEAT CAP AC 3 ; V OF VAPOR FREON- 1 T AT CONSTANT 
PRESSURE 

CP IS IN BTU/LBM DEB F . 

T IS IN DEB RANKING 






i 


FUNCTION CP ( T ) 
DT-1 . 

T3=T+DT 
P=PSAT ( T ) 

vo=vsm 

Vi=VG(T-5. > 
TOL-VO* 1 E--5 


2 V2»V0+<Vi~V0>/2. • 

IFt(Vl-VO) /2.LT.T0L) GO TO 1 
I F i FCPV3 • P , T3 , VO ) *FCPV3 < P , f 3 , V2 3 .. .8" . 0 . ) THEN 
VC"-V2 
ELSE 

V 1 ~-V2 |L®if 

SO TO 2 


Tiiiitt 


TS~T 

TOt =1 . E-4 

3 TP ::S TA+ (T8--TA) /2» 

JF C (TB-TA) /2, LT. TOL) BO TO 4 
in (mi T A ) — V3 ) * CVS ( TP > ~V3> . GT . O . 

, ELSE 

l» 

END IF 
GO TO 3 

4 TS«TP 
BB=.0O19 
C2=-33, 76999 
03=1,220367 
C5=~l . 478379E-4 
A 1=. 023815 

Cl ~2 . 798823E-4 
Dl =-2.1237341-7 
El =5. 999018E--1 1 
TC=848. 07 
AK*»-4.5 


3 THEN 
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* :V 


Mptllp vTS) -'«?B 
0 r re -*T3 


&• ! O — OK/"* C*TS '' O/TW’P 

fciali <C2/VB+, 5*C3/VB/VB+. 23*C5/V3**4 . ) * 

* 5 EXP (St ) * (ffil-l . ) -EXP (GTO) * . i I *. i^iis 



1 * < |T3^QSj^. 




; • • iCT ] UN R0PV3 WILL BE USED IN FUNCTION CP TO EVALUATE 
CP Zmm THE RELATIONSHIP BETWEEN CP AND CV 


FUNCTION FCPV3<P 5 T,V> .2 

R-. 073117 

BS=.0Ol9 

A2—3. 4-2*799 

£2**. 001318523 

C2=-35. 76999 

A3“-. 025341 

.87&121E-5 : >.i : T^ . : 

TT3=1. 220367 

A4- . 00 1 627277 

B4--1 . S05062E—6 v^pS.vi'C^ 

A5=-2. 35893E-5 
B5=2,44E303E-3 
. r;S«-: . 478379E-4 
A6~l . 057504E8 
E.6---9. 472I03E4 
TU*848.07 
AK=-4, 5 

•*P- E X J ‘ •, AK*T /7C > llli i ■ ' 

Vi..!= O— BB . 

F CPv3=R*T / VB+ < A2+B2*T+C2*PP ) /VB/VBf ( A3+B3*T+C3*PP> /VB**3. 
$ >:A4+34*T) /V3**4. 4 < A5+B5*T*C3*PP ) /VB**5. -P 
JF3V.8T. . i)CG TO 1 

FCF‘ - 3— FCf'v'3+ ! A6-;-B6*-T < /EXP (530. *V) 

1 RETURN 
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